notebook.community

Edit and run 





In [119]:



    


import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import linmix
from astropy.table import Table
import astropy.io.ascii as ascii
import corner
from lifelines import KaplanMeierFitter
from matplotlib import rc
from matplotlib.ticker import MultipleLocator

rc('font', **{'family': 'serif', 'serif': ['Computer Modern']})
rc('text', usetex=True)

#%matplotlib inline
%matplotlib qt5



    











In [229]:



    


Tab=Table.read("Table.fit")      #Taurus fluxes and disk mass Tab
Tab_Mass=Table.read("Table_Mass.fit")   #Taurus stellar masses (2-3 models) Tab
Lup_Tab=Table.read('Lupus_Tab.fit')      #Lupus Tab
ChaI_Tab=Table.read('ChaI_Tab.fit')
UppSco_Tab=Table.read('UppSco_Tab.fit')
#UppSco_Tab.show_in_browser()

#n=1     #Fmm
n=2     #Md

if n==2:
    T='20'   # T=20 K constant temperature
#    T='var'  # variable Temperature with luminosity

survey='Tau'   #Taurus
#survey='Lup'   #Lupus
#survey='ChaI'   #ChamaleonI
#survey='UppSco'   #Upper Scorpius OB

Introduction of variable delta --> if delta is True we have an observed data, if it is False we have an upper limit.





In [267]:



    


if survey=='Tau':
    delta=Tab['l']!='<'            #observed sources
if survey=='Lup':
    delta=Lup_Tab['l_Md']!='<'      #observed sources
if survey=='ChaI':
    delta=ChaI_Tab['l_Md']!='<'      #observed sources
if survey=='UppSco':
    delta=UppSco_Tab['l_Md']!='<'      #observed sources

notdelta = np.logical_not(delta)   #upperlimits



    











In [269]:



    


######################################
#        x and y data           #
if survey=='Tau':
    x=Tab['LogM*']
    x2=Tab_Mass['LogM*2']
    x3=Tab_Mass['LogM*3']   
    if n==1:
        y=np.log10(Tab['F1_3'])
    if n==2:
        if T=='var':
            y=Tab['LogMd']
        if T=='20':
            y=Tab['LogMd_20']    # Md in Eart Mass
        
if survey=='Lup':
    Md=Lup_Tab['Md']
    e_Md=Lup_Tab['err_Md']
    Md[notdelta]=3*e_Md[notdelta]    #Value for upperlimits   
    x=np.log10(Lup_Tab['M*'])
    y=np.log10(Md)  
    
if survey=='ChaI':
    Md=ChaI_Tab['Md']
    e_Md=ChaI_Tab['e_Md']   
    x=ChaI_Tab['logM*']
    y=np.log10(Md) 
    
if survey=='UppSco':
    Md=UppSco_Tab['Md_20']
    e_Md=UppSco_Tab['e_Md_20']    
    x=UppSco_Tab['logM*']
    y=np.log10(Md)  
    
############################################
#              Errors                #
if survey=='Tau': 
    ####errors on star mass (log scale)
    Dp_x=Tab['bM*_up']-Tab['LogM*']
    Dn_x=Tab['LogM*']-Tab['bM*_lo']
    Dp_x2=Tab_Mass['bM*_up2']-Tab_Mass['LogM*2']
    Dn_x2=Tab_Mass['LogM*2']-Tab_Mass['bM*_lo2']
    Dp_x3=Tab_Mass['bM*_up3']-Tab_Mass['LogM*3']
    Dn_x3=Tab_Mass['LogM*3']-Tab_Mass['bM*_lo3']
    ####errors on Continuum Flux (log scale)
    if n==1:
        Dp_y=np.log10(Tab['F1_3']+Tab['calibF'])-y
        Dn_y=y-np.log10(Tab['F1_3']-Tab['calibF'])
    ####errors on Disk mass (log scale)
    if n==2:
        if T=='var':
            Dp_y=Tab['Dp']   
            Dn_y=Tab['Dn']   
        if T=='20':
            Dp_y=Tab['Dp_20']   
            Dn_y=Tab['Dn_20'] 
        
if survey=='Lup':
    ecal_Md=np.sqrt((Lup_Tab['err_Md'])**2+(0.1*Lup_Tab['Md'])**2)    #quadratic sum of the calibration error on Md (10%)

    Dp_x=np.log10(Lup_Tab['M*']+Lup_Tab['err_M*'])-x  #superior error bar
    Dn_x=x-np.log10(Lup_Tab['M*']-Lup_Tab['err_M*'])
    Dp_y=np.log10(Md+ecal_Md)-y  #superior error bar
    Dn_y=y-np.log10(Md-ecal_Md)

if survey=='ChaI':
    ####errors on star mass (log scale)
    Dp_x=ChaI_Tab['up_logM*']-ChaI_Tab['logM*']
    Dn_x=ChaI_Tab['logM*']-ChaI_Tab['down_logM*']
    
    Dp_y=np.log10(Md+e_Md)-y  #superior error bar
    Dn_y=y-np.log10(Md-e_Md)   #inferior error bar
    
if survey=='UppSco':
    ####errors on star mass (log scale)
#    Dp_x=UppSco_Tab['up_logM*']-UppSco_Tab['logM*']
#    Dn_x=UppSco_Tab['logM*']-UppSco_Tab['down_logM*']
    Dp_x=UppSco_Tab['up_logM*']
    Dn_x=UppSco_Tab['down_logM*']
    
    Dp_y=np.log10(Md+e_Md)-y  #superior error bar
    Dn_y=y-np.log10(Md-e_Md)
    
#plt.errorbar(x,y,xerr=[Dn_x,Dp_x], yerr=[Dn_y,Dp_y], fmt='o', color='blue')    
#plt.show()



    











In [271]:



    


if survey=='Tau':
    x2sig=(Dp_x2+Dn_x2)/2.
    x3sig=(Dp_x3+Dn_x3)/2.
    
xsig=(Dp_x+Dn_x)/2.
ysig=(Dp_y+Dn_y)/2.

Run the Linmix code for the Bayesian linear regression analysis without considering the upper limits

lm = linmix.LinMix(x, y, xsig=xsig, ysig=ysig, K=2, nchains=30) lm.run_mcmc(miniter=5000, maxiter=100000, silent=True) #run the Markov Chain Monte Carlo
plt.errorbar(x,y,xerr=[Dn_x,Dp_x], yerr=[Dn_y,Dp_y], fmt='o', color='blue', alpha=0.5) #### for plot Fmm/Md vs M* for i in xrange(0, len(lm.chain), 150): xs = np.arange(-10,10) ys = lm.chain[i]['alpha'] + xs * lm.chain[i]['beta'] plt.plot(xs, ys, color='r', alpha=0.02) plt.xlim(-3,2) plt.ylim(-4,0) print lm.chain['alpha'].mean(), lm.chain['alpha'].std() print lm.chain['beta'].mean(), lm.chain['beta'].std() print lm.chain['sigsqr'].mean(), lm.chain['sigsqr'].std()

Considering the upper limits





In [124]:



    


lm_upp = linmix.LinMix(x, y, xsig=xsig, ysig=ysig, delta=delta, K=2, nchains=50)  #lin regression for Fmm/Md vs M*   
lm_upp.run_mcmc(miniter=5000, maxiter=100000, silent=True)    #run the Markov Chain Monte Carlo

Considering the mass derived from other models (BCAH98 and SDF00)





In [125]:



    


if survey=='Tau': 
    lm_upp2 = linmix.LinMix(x2, y, xsig=x2sig, ysig=ysig, delta=delta, K=2, nchains=50)    # lin regression for Fmm/Md vs (2)M*
    lm_upp2.run_mcmc(miniter=5000, maxiter=100000, silent=True)    #run the Markov Chain Monte Carlo



    











In [126]:



    


if survey=='Tau': 
    lm_upp3 = linmix.LinMix(x3, y, xsig=x3sig, ysig=ysig, delta=delta, K=2, nchains=50)    # lin regression for Fmm/Md vs (3)M*
    lm_upp3.run_mcmc(miniter=5000, maxiter=100000, silent=True)    #run the Markov Chain Monte Carlo



    











In [272]:



    


#### for plot Fmm/Md vs M*  ####
plt.errorbar(x[delta],y[delta] ,xerr=[Dn_x[delta],Dp_x[delta]], yerr=[Dn_y[delta],Dp_y[delta]], fmt='o', color='blue', alpha=0.5)
plt.errorbar(x[notdelta],y[notdelta],xerr=[Dn_x[notdelta],Dp_x[notdelta]],yerr=0.15,uplims=np.ones(sum(notdelta), dtype=bool), fmt='.', markersize='1',color='grey', alpha=0.3)

for i in xrange(0, len(lm_upp.chain), 80):
    xs = np.arange(-10,10)
    ys1 = lm_upp.chain[i]['alpha'] + xs * lm_upp.chain[i]['beta'] 
    plt.plot(xs,ys1, color='g', alpha=0.02)
    if survey=='Tau':
        plt.xlim(-3,2)
        plt.ylim(-1,3)
        if T=='var':
            plt.suptitle('Taurus', fontsize=20)
        if T=='20':
            plt.suptitle('Taurus T=20K', fontsize=20)
    if survey=='Lup':
        plt.xlim(-1.3,0.6)
        plt.ylim(-1.6,2.8)
        plt.suptitle('Lupus T=20K', fontsize=20)
    if survey=='ChaI':
        plt.xlim(-1.6,0.5)
        plt.ylim(-2,3)
        plt.suptitle('Chamaeleon I T=20K', fontsize=20)
    if survey=='UppSco':
        plt.xlim(-1,0.5)
        plt.ylim(-2,3)
        plt.suptitle('Upper Scorpius T=20K', fontsize=20)
    plt.xlabel('log(M$_{_*}$/M$_\odot$)', fontsize='16')
    plt.ylabel('log(M$_d$/M$_\odot$)', fontsize='16')
    
A=lm_upp.chain['alpha'].mean()
dA=lm_upp.chain['alpha'].std()    
B=lm_upp.chain['beta'].mean()
dB=lm_upp.chain['beta'].std()
sig=lm_upp.chain['sigsqr'].mean()
dsig=lm_upp.chain['sigsqr'].std()
corr=lm_upp.chain['corr'].mean()
dcorr=lm_upp.chain['corr'].std()
print A,dA
print B,dB
print np.sqrt(sig), dsig
print corr,dcorr



    


















    






1.41501990268 0.116091308345
1.71462305971 0.206192397981
0.693313078253 0.0785725160201
0.677971179724 0.0602435117904
208

















In [273]:



    


if survey=='Tau': 
    #### for plot Fmm/Md vs M*(2)  ####

    plt.errorbar(x2[delta],y[delta] ,xerr=[Dn_x2[delta],Dp_x2[delta]], yerr=[Dn_y[delta],Dp_y[delta]], fmt='o', color='blue', alpha=0.5)
    plt.errorbar(x2[notdelta],y[notdelta],xerr=[Dn_x2[notdelta],Dp_x2[notdelta]],yerr=0.15,uplims=np.ones(sum(notdelta), dtype=bool), fmt='.', markersize='1',color='grey', alpha=0.3)

    for i in xrange(0, len(lm_upp.chain), 80):
        xs = np.arange(-10,10)
        ys2 = lm_upp2.chain[i]['alpha'] + xs * lm_upp2.chain[i]['beta'] 
        plt.plot(xs, ys2, color='g', alpha=0.02)
        plt.xlim(-3,2)
        plt.ylim(-4,0)

    A2=lm_upp2.chain['alpha'].mean()
    dA2=lm_upp2.chain['alpha'].std()    
    B2=lm_upp2.chain['beta'].mean()
    dB2=lm_upp2.chain['beta'].std()
    sig2=lm_upp2.chain['sigsqr'].mean()
    dsig2=lm_upp2.chain['sigsqr'].std()
    corr2=lm_upp2.chain['corr'].mean()
    dcorr2=lm_upp2.chain['corr'].std()
    print A2, dA2
    print B2, dB2
    print np.sqrt(sig2), dsig2
    print corr2, dcorr2
    
    #### for plot Fmm/ Md vs M*(3) ####
    plt.errorbar(x3[delta],y[delta] ,xerr=[Dn_x3[delta],Dp_x3[delta]], yerr=[Dn_y[delta],Dp_y[delta]], fmt='o', color='blue', alpha=0.5)
    plt.errorbar(x3[notdelta],y[notdelta],xerr=[Dn_x3[notdelta],Dp_x3[notdelta]],yerr=0.15,uplims=np.ones(sum(notdelta), dtype=bool), fmt='.', markersize='1',color='grey', alpha=0.3)

    for i in xrange(0, len(lm_upp.chain), 50):
        xs = np.arange(-10,10)
        ys3 = lm_upp3.chain[i]['alpha'] + xs * lm_upp3.chain[i]['beta'] 
        plt.plot(xs, ys3, color='g', alpha=0.02)
        plt.xlim(-3,2)
        plt.ylim(-4,0)

    A3=lm_upp3.chain['alpha'].mean()
    dA3=lm_upp3.chain['alpha'].std()    
    B3=lm_upp3.chain['beta'].mean()
    dB3=lm_upp3.chain['beta'].std()
    sig3=lm_upp3.chain['sigsqr'].mean()
    dsig3=lm_upp3.chain['sigsqr'].std()
    corr3=lm_upp3.chain['corr'].mean()
    dcorr3=lm_upp3.chain['corr'].std()

    print A3, dA3
    print B3, dB3
    print np.sqrt(sig3), dsig3
    print corr3, dcorr3



    


















    






0.938309643819 0.0700068750686
1.27862007608 0.143457520857
0.693316405819 0.0764394224461
0.665678079897 0.0472017924567
1.18200638608 0.0819571716715
1.55417410693 0.162467314231
0.666809967971 0.0698483149419
0.692224713257 0.0486874554815
208

















In [274]:



    


if survey=='Tau':
    fig, ax = plt.subplots(1, 3, figsize=(5, 5) ,sharey=True)
    plt.subplots_adjust(wspace=0.1,hspace=0.01)
    Ms=(x,x2,x3)
    Dp_Ms=(Dp_x,Dp_x2,Dp_x3)
    Dn_Ms=(Dn_x,Dn_x2,Dn_x3)
    model=('DM97','BCAH98','SDF00')
    color=('red','green','cyan')

    ax[0].set_xlabel('log(M$_{_*}$/M$_\odot$)', fontsize='16')

    for j in range(3):
        if n==1:
            ax[j].set_ylim(-3.5,0.1)
            ax[j].text(-2.2, -0.2, model[j],color=color[j], fontsize=15)
        if n==2:
            ax[j].set_ylim(-4.1,-0.5)  
            ax[j].text(-2.2, -0.8, model[j],color=color[j], fontsize=15)
        ax[j].set_xlim(-2.5,1.5)
        ax[j].errorbar(Ms[j][delta],y[delta] ,xerr=[Dn_Ms[j][delta],Dp_Ms[j][delta]], yerr=[Dn_y[delta],Dp_y[delta]], fmt='o', color='blue', alpha=0.5)
        ax[j].errorbar(Ms[j][notdelta],y[notdelta],xerr=[Dn_Ms[j][notdelta],Dp_Ms[j][notdelta]],yerr=0.15,uplims=np.ones(sum(notdelta), dtype=bool), fmt='.', markersize='1',color='grey', alpha=0.3)
        ax[j].tick_params(axis='x', labelsize=12)
        ax[j].tick_params(axis='y', labelsize=12)
    
    for i in xrange(0, len(lm_upp.chain), 50):
        xs = np.arange(-10,10)
        ys = lm_upp.chain[i]['alpha'] + xs * lm_upp.chain[i]['beta'] 
        ys2 = lm_upp2.chain[i]['alpha'] + xs * lm_upp2.chain[i]['beta'] 
        ys3 = lm_upp3.chain[i]['alpha'] + xs * lm_upp3.chain[i]['beta'] 
        ax[0].plot(xs, ys, color='pink', alpha=0.02)
        ax[1].plot(xs, ys2, color='g', alpha=0.02)
        ax[2].plot(xs, ys3, color='cyan', alpha=0.02)

    if n==1:
        ax[0].set_ylabel('log(F$_{mm}$/Jy)', fontsize='16')
        plt.suptitle('F$_{mm}$ vs M$_{_*}$', fontsize=20)
    if n==2:
        ax[0].set_ylabel('log(M$_d$/M$_\odot$)', fontsize='16')
        plt.suptitle('M$_d$ vs M$_{_*}$', fontsize=20)



    


















    






208

writing output of regression analysis in an ascii file





In [221]:



    


if survey=='Tau':
    if n==1:
        ascii.write(lm_upp.chain[['alpha', 'beta', 'sigsqr','mu0', 'usqr', 'wsqr','ximean', 'xisig', 'corr']],
                    'Fmm_vs_Ms.pyout', overwrite=True)
        ascii.write(lm_upp2.chain[['alpha', 'beta', 'sigsqr','mu0', 'usqr', 'wsqr','ximean', 'xisig', 'corr']],
                    'Fmm_vs_Ms2.pyout', overwrite=True)
        ascii.write(lm_upp3.chain[['alpha', 'beta', 'sigsqr','mu0', 'usqr', 'wsqr','ximean', 'xisig', 'corr']],
                    'Fmm_vs_Ms3.pyout', overwrite=True)
    if n==2:
        if T=='var':
            ascii.write(lm_upp.chain[['alpha', 'beta', 'sigsqr','mu0', 'usqr', 'wsqr','ximean', 'xisig', 'corr']],
                    'Md_vs_Ms.pyout', overwrite=True)
            ascii.write(lm_upp2.chain[['alpha', 'beta', 'sigsqr','mu0', 'usqr', 'wsqr','ximean', 'xisig', 'corr']],
                    'Md_vs_Ms2.pyout', overwrite=True)
            ascii.write(lm_upp3.chain[['alpha', 'beta', 'sigsqr','mu0', 'usqr', 'wsqr','ximean', 'xisig', 'corr']],
                    'Md_vs_Ms3.pyout', overwrite=True)
        if T=='20':
            ascii.write(lm_upp.chain[['alpha', 'beta', 'sigsqr','mu0', 'usqr', 'wsqr','ximean', 'xisig', 'corr']],
                    'Md_20_vs_Ms.pyout', overwrite=True)
            ascii.write(lm_upp2.chain[['alpha', 'beta', 'sigsqr','mu0', 'usqr', 'wsqr','ximean', 'xisig', 'corr']],
                    'Md_20_vs_Ms2.pyout', overwrite=True)
            ascii.write(lm_upp3.chain[['alpha', 'beta', 'sigsqr','mu0', 'usqr', 'wsqr','ximean', 'xisig', 'corr']],
                    'Md_20_vs_Ms3.pyout', overwrite=True)
if survey=='Lup':
    ascii.write(lm_upp.chain[['alpha', 'beta', 'sigsqr','mu0', 'usqr', 'wsqr','ximean', 'xisig', 'corr']],
                    'Lup_Md_vs_Ms.pyout', overwrite=True)
if survey=='ChaI':
    ascii.write(lm_upp.chain[['alpha', 'beta', 'sigsqr','mu0', 'usqr', 'wsqr','ximean', 'xisig', 'corr']],
                    'ChaI_Md_vs_Ms.pyout', overwrite=True)
if survey=='UppSco':
    ascii.write(lm_upp.chain[['alpha', 'beta', 'sigsqr','mu0', 'usqr', 'wsqr','ximean', 'xisig', 'corr']],
                    'UppSco_Md_vs_Ms.pyout', overwrite=True)

Corner plots

if n==1: pyout = ascii.read('Fmm_vs_Ms.pyout') if n==2: pyout = ascii.read('Md_vs_Ms.pyout') fig, axarr = plt.subplots(4, 4, figsize=(10, 10)) corner.corner(np.array([pyout['alpha'], pyout['beta'], pyout['sigsqr'],pyout['corr'] ]).T, labels=[r"$\alpha$", r"$\beta$", r"$\sigma^2$",r"corr"], range=[0.99]*4, plot_datapoints=False,show_titles=True, fig=fig) fig.subplots_adjust(bottom=0.065, left=0.07) plt.show()

Posterior Distributions





In [276]:



    


if survey=='Tau':
   
    if n==1:
        pyout1 = ascii.read('Fmm_vs_Ms.pyout')
        pyout2 = ascii.read('Fmm_vs_Ms2.pyout')
        pyout3 = ascii.read('Fmm_vs_Ms3.pyout')
    if n==2:
        if T=='var':
            pyout1 = ascii.read('Md_vs_Ms.pyout')
            pyout2 = ascii.read('Md_vs_Ms2.pyout')
            pyout3 = ascii.read('Md_vs_Ms3.pyout')
        if T=='20':
            pyout1 = ascii.read('Md_20_vs_Ms.pyout')
            pyout2 = ascii.read('Md_20_vs_Ms2.pyout')
            pyout3 = ascii.read('Md_20_vs_Ms3.pyout')

    a1=pyout1['alpha']
    b1=pyout1['beta']
    s1=np.sqrt(pyout1['sigsqr'])
    corr1=pyout1['corr']

    a2=pyout2['alpha']
    b2=pyout2['beta']
    s2=np.sqrt(pyout2['sigsqr'])
    corr2=pyout2['corr']

    a3=pyout3['alpha']
    b3=pyout3['beta']
    s3=np.sqrt(pyout3['sigsqr'])
    corr3=pyout3['corr']

if survey=='Lup':
    pyout1 = ascii.read('Lup_Md_vs_Ms.pyout')
    a1=pyout1['alpha']
    b1=pyout1['beta']
    s1=np.sqrt(pyout1['sigsqr'])
    corr1=pyout1['corr']

if survey=='ChaI':
    pyout1 = ascii.read('ChaI_Md_vs_Ms.pyout')
    a1=pyout1['alpha']
    b1=pyout1['beta']
    s1=np.sqrt(pyout1['sigsqr'])
    corr1=pyout1['corr']

if survey=='UppSco':
    pyout1 = ascii.read('UppSco_Md_vs_Ms.pyout')
    a1=pyout1['alpha']
    b1=pyout1['beta']
    s1=np.sqrt(pyout1['sigsqr'])
    corr1=pyout1['corr']



    











In [278]:



    


w1 = np.ones_like(a1)/float(len(a1))
fig, ax = plt.subplots(1, 4, figsize=(5, 5) ,sharey=True) plt.subplots_adjust(wspace=0,hspace=0.01) label=('alpha','beta','RMS','corr') for j in range(4): ax[j].tick_params(axis='x', labelsize=12) ax[j].tick_params(axis='y', labelsize=12) ax[j].set_ylim(0,9) ax[j].set_xlabel(label[j],fontsize=12) ax[0].hist(a1,20,normed=True, edgecolor='red', hatch="///",histtype='step') ax[1].hist(b1,20,normed=True,edgecolor='red', hatch="///",histtype='step') ax[2].hist(s1,20,normed=True,edgecolor='red', hatch="///",histtype='step') ax[3].hist(corr1,20,normed=True,edgecolor='red', hatch="///",histtype='step') plt.show() 




In [280]:



    


fig1, ax1 = plt.subplots(1, 4, figsize=(5, 5) ,sharey=True)
plt.subplots_adjust(wspace=0,hspace=0.01)
label=('alpha','beta','RMS','corr')

for j in range(4):
    ax1[j].tick_params(axis='x', labelsize=12)
    ax1[j].tick_params(axis='y', labelsize=12)
    ax1[j].set_ylim(0,0.45)
    ax1[j].set_xlabel(label[j],fontsize=12)
    
if survey=='Tau':
    if n==1:
        plt.suptitle('Taurus F$_{mm}$',fontsize=20)
        ax1[0].set_ylabel('P(X|M$_{_*}$,F$_{mm}$)', fontsize=16)
        ax1[0].set_xlim(-2.6,-0.6)
        ax1[1].set_xlim(0.4,2.9)
        ax1[2].set_xlim(0.3,1.05)
        ax1[3].set_xlim(-0.1,1.05)
        bins1=np.arange(-2.6,-0.6,0.05)
        bins2=np.arange(0.4,2.9,0.05)
    if n==2:
        if T=='var':
            plt.suptitle('Taurus M$_d$', fontsize=20)    
            ax1[0].set_ylabel('P(X|M$_{_*}$,M$_d$)', fontsize=16)
            ax1[0].set_xlim(-3.6,-1.2)
            ax1[1].set_xlim(-0.2,2.8)
            ax1[2].set_xlim(0.3,1.05)
            ax1[3].set_xlim(-0.1,1.05)    
            bins1=np.arange(-3.6,-1.2,0.05)
            bins2=np.arange(-0.2,2.8,0.05)
        if T=='20':
            plt.suptitle('Taurus M$_d$ 20K', fontsize=20)    
            ax1[0].set_ylabel('P(X|M$_{_*}$,M$_d$)', fontsize=16)
            ax1[0].set_xlim(0.5,2.4)
            ax1[1].set_xlim(0.5,3.1)
            ax1[2].set_xlim(0.3,1.1)
            ax1[3].set_xlim(-0.1,1.05)       
            bins1=np.arange(0.5,2.4,0.05)
            bins2=np.arange(0.5,3.1,0.05)   
    
    bins3=np.arange(0.3,1.05,0.05)
    bins4=np.arange(-0.1,1.05,0.05)
    
if survey=='Lup':
    plt.suptitle('Lupus M$_d$', fontsize=20)    
    ax1[0].set_ylabel('P(X|M$_{_*}$,M$_d$)', fontsize=16)
    ax1[0].set_xlim(0.5,2.4)
    ax1[1].set_xlim(0.5,3.1)
    ax1[2].set_xlim(0.3,1.1)
    ax1[3].set_xlim(-0.1,1.05)    
    bins1=np.arange(0.5,2.4,0.05)
    bins2=np.arange(0.5,3.1,0.05)    
    bins3=np.arange(0.3,1.1,0.05)
    bins4=np.arange(-0.1,1.05,0.05)

if survey=='ChaI':
    plt.suptitle('ChameleonI M$_d$', fontsize=20)    
    ax1[0].set_ylabel('P(X|M$_{_*}$,M$_d$)', fontsize=16)
    ax1[0].set_xlim(0.2,2)
    ax1[1].set_xlim(0.7,3.2)
    ax1[2].set_xlim(0.3,1.1)
    ax1[3].set_xlim(-0.1,1.05)    
    bins1=np.arange(0.2,2,0.05)
    bins2=np.arange(0.7,3.2,0.05)    
    bins3=np.arange(0.3,1.1,0.05)
    bins4=np.arange(-0.1,1.05,0.05)
    
if survey=='UppSco':
    plt.suptitle('Upper Scorpius OB M$_d$', fontsize=20)    
    ax1[0].set_ylabel('P(X|M$_{_*}$,M$_d$)', fontsize=16)
    ax1[0].set_xlim(0,1.7)
    ax1[1].set_xlim(1,3.8)
    ax1[2].set_xlim(0.3,1.1)
    ax1[3].set_xlim(-0.1,1.05)    
    bins1=np.arange(0,1.7,0.05)
    bins2=np.arange(1,3.8,0.05)    
    bins3=np.arange(0.3,1.1,0.05)
    bins4=np.arange(-0.1,1.05,0.05)
    
    
ax1[0].hist(a1,bins=bins1 ,weights=w1, edgecolor='red', hatch="\\\ ",histtype='step',label='DM97')
ax1[1].hist(b1,bins=bins2 ,weights=w1,edgecolor='red', hatch="\\\ ",histtype='step',label='DM97')
ax1[2].hist(s1,bins=bins3 ,weights=w1,edgecolor='red', hatch="\\\ ",histtype='step',label='DM97')
ax1[3].hist(corr1,bins=bins4 ,weights=w1,edgecolor='red', hatch="\\\ ",histtype='step',label='DM97')

if survey=='Tau':
    ax1[0].hist(a2,bins=bins1 ,weights=w1,edgecolor='green', hatch="//",histtype='step',label='BCAH98')
    ax1[0].hist(a3,bins=bins1 ,weights=w1,edgecolor='cyan', hatch="\\ ",histtype='step',label='SDF00')
    ax1[1].hist(b2,bins=bins2 ,weights=w1,edgecolor='green', hatch="//",histtype='step',label='BCAH98')
    ax1[1].hist(b3,bins=bins2 ,weights=w1,edgecolor='cyan', hatch="\\ ",histtype='step',label='SDF00')
    ax1[2].hist(s2,bins=bins3 ,weights=w1,edgecolor='green', hatch="//",histtype='step',label='BCAH98')
    ax1[2].hist(s3,bins=bins3 ,weights=w1,edgecolor='cyan', hatch='\\ ',histtype='step',label='SDF00')
    ax1[3].hist(corr2,bins=bins4 ,weights=w1,edgecolor='green', hatch="//",histtype='step',label='BCAH98')
    ax1[3].hist(corr3,bins=bins4 ,weights=w1,edgecolor='cyan', hatch="\\ ",histtype='step',label='SDF00')
    plt.legend()

plt.show()

Cumulative distribution from the Kaplan Meier estimator





In [296]:



    


#plot='ratios'
plot='diskmass'

if plot=='ratios':           #survival distrib for ratios Mdisk/Mstellar_host
    if survey=='Tau':
        if T=='var':
            data1 = (10**Tab['LogMd'])/(10**Tab['LogM*'])    
            data2 = (10**Tab['LogMd'])/(10**Tab_Mass['LogM*2'])    
            data3 = (10**Tab['LogMd'])/(10**Tab_Mass['LogM*3'])  
        if T=='20':
            data1 = (10**Tab['LogMd_20'])/(10**Tab['LogM*'])    
            data2 = (10**Tab['LogMd_20'])/(10**Tab_Mass['LogM*2'])    
            data3 = (10**Tab['LogMd_20'])/(10**Tab_Mass['LogM*3'])  
    if survey=='Lup':
        data1 = Md/(Lup_Tab['M*'])  
    if survey=='ChaI':
        data1 = Md/(10**ChaI_Tab['logM*'])          
    if survey=='UppSco':
        data1 = Md/(10**UppSco_Tab['logM*'])  
        
if plot=='diskmass':        #survival distrib for Disk mass
    if survey=='Tau':
        if T=='var':
            data1 = 10**Tab['LogMd'] 
            data2 = 10**Tab['LogMd']   
            data3 = 10**Tab['LogMd']
        if T=='20':
            data1 = 10**Tab['LogMd_20'] 
            data2 = 10**Tab['LogMd_20']   
            data3 = 10**Tab['LogMd_20']
    if survey=='Lup':
        data1 = Md*1. 
    if survey=='ChaI':
        data1 = Md*1. 
    if survey=='UppSco':
        data1 = Md*1. 

C=delta
kmf_1= KaplanMeierFitter()
kmf_1.fit(data1, event_observed=C,left_censorship=True,label="kmf",alpha=0.5)
s1=1-kmf_1.cumulative_density_.values.ravel()
t1=kmf_1.timeline
Dp1=kmf_1.confidence_interval_['kmf_upper_0.50']-(1-s1)
Dn1=(1-s1)-kmf_1.confidence_interval_['kmf_lower_0.50']
ax = plt.subplot(111)

if survey=='Tau':
    kmf_2= KaplanMeierFitter()
    kmf_3= KaplanMeierFitter()
    
    kmf_2.fit(data2, event_observed=C,left_censorship=True,label="kmf",alpha=0.5)
    kmf_3.fit(data3, event_observed=C,left_censorship=True,label="kmf",alpha=0.5)

    s2=1-kmf_2.cumulative_density_.values.ravel()
    s3=1-kmf_3.cumulative_density_.values.ravel()
    t2=kmf_2.timeline
    t3=kmf_3.timeline

    Dp2=kmf_2.confidence_interval_['kmf_upper_0.50']-(1-s2)
    Dp3=kmf_3.confidence_interval_['kmf_upper_0.50']-(1-s3)
    Dn2=(1-s2)-kmf_2.confidence_interval_['kmf_lower_0.50']
    Dn3=(1-s3)-kmf_3.confidence_interval_['kmf_lower_0.50']

    Dp=(Dp1,Dp2,Dp3)
    Dn=(Dn1,Dn2,Dn3)
    model=('DM97','BCAH98','SDF00')
    color=('red','green','blue')
    t=(t1,t2,t3)
    s=(s1,s2,s3)

    for i in range(3):
        ax.errorbar(t[i],s[i],yerr=[Dn[i],Dp[i]],drawstyle='steps-post',linewidth=3,
                color=color[i],label=model[i],elinewidth=1)

if survey=='Tau':
    label='Taurus'
    plt.xlim(5e-1, 3e2)
if survey=='Lup':
    label='Lupus'
    plt.xlim(5e-2, 3e2)
if survey=='ChaI':
    label='Chamaeleon I'
    plt.xlim(5e-2, 3e2)
if survey=='UppSco':
    label='Upper Scorpius'
    plt.xlim(5e-2, 3e2)
    
ax.errorbar(t1,s1,yerr=[Dn1,Dp1],drawstyle='steps-post',linewidth=3,
                color='red',label=label,elinewidth=1)

minorLocator   = MultipleLocator(0.05)
ax.yaxis.set_minor_locator(minorLocator)

plt.legend()
plt.xscale('log')
plt.ylim(0, 1)
if plot=='ratios':
    plt.ylabel('f($>$M$_d$/M$_{_*}$)', fontsize=16)
    plt.xlabel('M$_d$/M$_{_*}$',fontsize=16)
    
if plot=='diskmass':
    plt.ylabel('f($>$M$_d$)', fontsize=16)
    plt.xlabel('M$_d$',fontsize=16)

#print "Median survival1 :", kmf_1.median_
# print "Median survival2 :", kmf_2.median_
# print "Median survival3 :", kmf_3.median_
plt.show()

Cumulative Distributions of the disk to mass ratios





In [146]:



    


sorted_data1=np.sort(data1)
p1=np.arange(len(sorted_data1))/float(len(sorted_data1)-1)
plt.step(sorted_data1,1-p1, 'red')    #1-p is the survival function curve

if survey=='Tau':
    sorted_data2=np.sort(data2)
    sorted_data3=np.sort(data3)
    p2=np.arange(len(sorted_data2))/float(len(sorted_data2)-1)
    p3=np.arange(len(sorted_data3))/float(len(sorted_data3)-1)
    plt.step(sorted_data2,1-p2,'green')    #1-p is the survival function curve
    plt.step(sorted_data3,1-p3, 'blue')    #1-p is the survival function curve

plt.xscale('log')
plt.ylabel('f($>$M$_d$/M$_{_*}$)', fontsize=16)
plt.xlabel('M$_d$/M$_{_*}$',fontsize=16)
plt.show()

PROVA





In [148]:



    


ax = plt.subplot(111)
Tab4=Table.read("Tab4_totbin.fit")
Tab4.remove_rows([8,10,31,34,76,94,103,110,118,125,126,174,180,195,204,220]) 

Fmm=Tab['F1_3']
SpT=Tab['SpT']
Ts=Tab4['logT_']
'O'<'A'<'F'<'G'<'K'<'M'
ran=np.arange(3.5241,3.6238,0.0001)

d= SpT<'K6'
d1=SpT>'M6'
d2=Ts>3.6238   #< K6
d3=Ts<3.4757    #> M6
d4=np.empty(len(Ts),dtype=bool)
d5=np.empty(len(Ts),dtype=bool)
for i in range(0,len(Ts)):
    d4[i]= 3.5241<Ts[i]<3.6238
    d5[i]= 3.4757<Ts[i]<3.5241
kmf1 = KaplanMeierFitter()
kmf2 = KaplanMeierFitter()
kmf3 = KaplanMeierFitter()
    
C=delta
# kmf1.fit(Fmm[d], event_observed=C[d],alpha=0.1,label='$<$ K6').plot(ax=ax,color='purple')
# kmf2.fit(Fmm[d1], event_observed=C[d1],alpha=0.1,label='$>$ M6').plot(ax=ax,color='red')
kmf1.fit(Fmm[d2], event_observed=C[d2],alpha=0.1,label='$<$ K6(T)').plot(ax=ax,color='purple')
kmf2.fit(Fmm[d3], event_observed=C[d3],alpha=0.1,label='$>$ M6(T)').plot(ax=ax,color='red')
kmf2.fit(Fmm[d4], event_observed=C[d4],alpha=0.1,label=' K6(T)-M3.5(T)').plot(ax=ax,color='blue')
kmf2.fit(Fmm[d5], event_observed=C[d5],alpha=0.1,label=' M3.5(T)-M6(T)').plot(ax=ax,color='orange')
plt.xscale('log')
plt.ylim(0,1)
plt.xlim(2e-4, 1)
plt.ylabel('f($>F_{mm}$)', fontsize=16)
plt.xlabel('F$_{mm}$',fontsize=16)
major_ticks = np.arange(0, 1, 0.1)                                              
minor_ticks = np.arange(0, 1, 0.05) 
ax.set_yticks(minor_ticks, minor=True)  
# ax.tick_params(axis='y',right='on')
print SpT[11],Ts[11]



    


















    






M3.5 3.5241










    






/home/mdesimon/bin/anaconda2/lib/python2.7/site-packages/matplotlib/cbook.py:2649: UserWarning: Saw kwargs [u'c', u'color'] which are all aliases for u'color'.  Kept value from u'color'
  seen=seen, canon=canonical, used=seen[-1]))

















In [150]:



    


plt.plot(Ts[delta],Fmm[delta],'.r')
plt.plot(Ts[notdelta],Fmm[notdelta],'.g',alpha=0.5)
plt.yscale('log')
plt.ylim(4e-4,1)



    


















    
Out[150]:








(0.0004, 1)

















In [ ]:

Content source: martadesimone/Protoplanetarydisks

Similar notebooks:

notebook.community | gallery | about